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Abstract 

We study the anomalous flux ratio which is observed in some four-image lens sys- 
tems, where the source lies close to a fold caustic. In this case two of the images 
are close to the critical curve and their flux ratio should be equal to unity, instead 
in several cases the observed value differs significantly. The most plausible solution 
is to invoke the presence of substructures, as for instance predicted by the Cold 
Dark Matter scenario, located near the two images. In particular, we analyze the 
two fold lens systems PG1115+080 and B1555+375, for which there are not yet 
satisfactory models which explain the observed anomalous flux ratios. We add to 
a smooth lens model, which reproduces well the positions of the images but not 
the anomalous fluxes, one or two substructures described as singular isothermal 
spheres. For PG1115+080 we consider a smooth model with the influence of the 
group of galaxies described by a SIS and a substructure with mass ~ 10 5 M@ as 
well as a smooth model with an external shear and one substructure with mass 
~ 10 8 Mq. For B1555+375 either a strong external shear or two substructures with 
mass ~ 10 7 M Q reproduce the data quite well. 

Key words: cosmology: theory - dark matter - gravitational lensing - galaxies: 
haloes - substructures 



1 Introduction 



The standard lens models, although reproduce in general the relative positions 
of the images quite accurately, often have difficulties explaining the relative 
fluxes of multiply-imaged sources. Several possible explanations have been 
considered in the literature, the most plausible being that the lensing potential 
of real galaxies are not fully described by the simple lens models used to 
compute the lens characteristics. The most often invoked solution is to consider 
additional small-scale structures, which if located near the images can modify 
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significantly the observed flux ratio between different images, in particular the 
so-called cusp or fold relations. 

The presence of substructures is naturally expected within the Cold Dark 
Matter (CDM) model, which has been successful in explaining a large variety 
of observational results like the large scale structure of galaxies on scales larger 
than 1 Mpc or the fluctuations of the CMB ( Spergel et al 2003). However, 
one of the predictions of this scenario is a distribution of matter, with a large 
number of small-mass compact dark matter (DM subhalos) halos, both within 
virialized regions of larger halos (Moore et al. 1999, Klypin et al. 1999) and 
in the field (DM extragalactic halos, Metcalf 2005). At the same time, the 
observed number of dwarf galaxy satellites in the Local Group is more than 
an order of magnitude smaller than expected. Many theoretical studies suggest 
models to reduce the abundance of substructure or to suppress star formation 
in small clumps via astrophysical mechanisms (as feedback), making them 
dark (Bullock, Kravtosov & Weinberg 2000, Kravtsov et al. 2004, Moore et 
al. 2006). Anyway, if the CDM paradigm is correct, we expect ~ 10 — 15% of 
the mass of a present-day galactic halo (~ 1O 12 M ) within the virial radius to 
be in substructures with mass > 10 7 M Q . Thus in the CDM model anomalous 
flux ratios should be common. 

At present, the only way to detect these subclumps is through gravitational 
lensing, which is directly sensitive to the mass. This is because substructures 
(like globular clusters, gas clouds or satellite galaxies) can strongly modify 
the fluxes of lensed images relative to those predicted by smooth lens models. 
Even a clump as small as a star can perturb the image of small sources (~ 100 
AU), but in this case we would see a microlensing effect such as variability in 
the image brightness with a time scale of order months. Some authors tried 
to fit systems (for example the radio system B1422+231) with anomalous flux 
ratios using smooth lens and multipole models: although it seems necessary to 
investigate further whether the multipoles can fit the lens configurations, these 
methods are not exhaustive (Evans & Witt 2003, Kochanek & Dalai 2004, 
Congdon & Keeton 2005) and it is not yet possible to conclude that the mul- 
tipole approach can explain the anomalous flux ratios. Indeed, Kochanek & 
Dalai 2004 showed that the flux anomaly distributions display the characteris- 
tic demagnifications of the brightest saddle point relative to the other images 
expected for low optical depth substructure, which they conclude cannot be 
mimicked by problems in the "macro" models for the gravitational potential 
of the lens galaxy. Mao & Schneider (1998), Keeton (2001), Metcalf & Madau 
(2001), Bradac et al. (2002), Dobler &Keeton (2006) noted that a simple way 
of solving the puzzle was to put a satellite near the images, and they found 
that this could explain the anomaly in B1422+231. 

Generally the flux ratios between the images do not depend on wavelength 
since they are independent of the intrinsic flux and variability of the source 
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(Keeton et al. 1997, Mao & Schneider 1998, Keeton 2001, Metcalf & Zhao 
2002). Such discrepancies would probably be due to sub-lensing or microlens- 
ing effects. On the other hand when modeling a multiple QSO lens system, 
one can either include or disregard the flux ratios of the images. As pointed 
out by Chang & Refsdal (1979) ( and in the following by several authors, 
as for instance Metcalf 2005, Keeton et al. 2005, Mortonson et al. 2005) we 
have to pay attention to the fact that the projected (on the lens plane) sizes 
of the optical continuum emitting regions of QSOs are expected to be of the 
same order as the Einstein radius of a star in the lens galaxy (~ 100AU), so 
that the optical magnitudes may well be affected by gravitational microlensing 
(see Metcalf et al. 2004, Mortonson et al. 2005, and references therein), even 
if averaged over long periods of time. The radio and mid-IR regions, when pro- 
jected on the lens plane, are typically of the order of 10 pc and change in their 
magnification should be dominated by larger scales than stars (Metcalf 2005, 
Chiba et al. 2005). If the lens galaxies contain substructures with an Einstein 
radius comparable or greater than the projected size of the radio component 
(corresponding for the substructure to masses > 1O 8 M ) we should see image 
splitting and distortions if they lie close enough to the images (see section 
3.2, Wambsganss & Paczyhski 1994, Mao & Schneider 1998), which have not 
yet been detected (this might be the case for B0128+437, Biggs et al. 2004). 

The existence of anomalous fluxes in many lens systems has been known since 
some time. The substructure lensing effects have been studied by considering 
single substructures (Mao & Schneider 1998, Metcalf & Madau 2001), by as- 
suming a statistically distributed sample of substructures (Chiba 2002, Chen 
et al. 2003, Keeton et al. 2005) or by simulations (Amara et al. 2006, Maccio 
et al. 2005)13 An explanation by lensing of substructures for the anomalous 
flux ratio for the PG1115+080 system has already been considered (Chiba 
2002, Dalai & Kochanek 2002, Chen et al. 2003, Kochanek & Dalai 2004, 
Keeton, Gaudi, Petters 2005), however, mainly within a statistical treatment 
to determine whether a plausible collection of mass clumps could explain the 
strange flux ratio. As pointed out by Chiba (2002), even if it seems difficult to 
reproduce anomalous flux ratios with CDM subhalos, he concluded that the 
main role in reproducing the observed flux ratio is played by the one satellite 
which is located in the vicinity of an image (either Al or A2 in PG1115+80). 

In this paper we analyze in detail the lens system by adding one or two sub- 
clumps nearby one of the images and by solving the lens equation. In section 
2 we review the main observation and analysis done so far on the two lens sys- 
tems PG1115+080 and B1555+375. In section 3 we briefly recall the relevant 
formalism for gravitational lensing and how we proceed when we consider a 
lens model with a perturbation induced by one or more substructures. Assum- 
ing a SIS model for the substructures we can then get an estimate on their 



These works deal with violations of the cusp relation. 
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position and Einstein radius (or mass) such as to modify the flux of the image 
pair near the critical curve due to a source located close to a fold. In section 
4 we present the numerical simulations and fits to the two considered lens 
systems. We conclude with a short summary and discussion of our results in 
section 5. 



2 About PG1115+080 and B1555+375 

PG1115+080 is the second gravitationally lensed quasar which was discovered 
(Weymann et al. 1980, Impey et al. 1998 and references therein). The source 
is at redshift z s = 1.722 and the lens galaxy at z\ = 0.310. It is an optically 
selected, radio-quiet quasar. Hege et al. (1981) first resolved the four quasar 
images (a close pair A2/A1, B and C), confirming the early model of Young 
et al. (1981) that the lens is a five- images system, one image being hidden by 
the core of the lens galaxy. 

Young et al. (1981) noted that the lens galaxy seems to be part of a small 
group centered to the southwest of the lens, with a velocity dispersion of 
approximately 270 ± 70 km s _1 based only on four galaxy redshifts ( Tonry 
1998 also confirmed the group velocity dispersion by using five galaxies and 
getting a value of 326 kms^ 1 ). The group is an essential component of any 
model to successfully fit the lens constraints (Keeton et al. 1997; Schechter 
et al. 1997). Also two time delays between the images were determined by 
Schechter et al. (1997) and confirmed by Barkana (1997). Their results were 
analyzed by Keeton & Kochanek (1997) and Courbin et al. (1997) by assuming 
power law mass distributions and leading to a value of H = 53~!l\ 5 km s _1 
Mpc -1 , with comparable contributions to the uncertainties both from the time 
delay measurements and the models. Recently, Read et al. 2007, assuming a 
non parametric mass distribution, found H = 64jJ km s -1 Mpc -1 consistent 
with the currently accepted value of about 70 km s -1 Mpc -1 . 

We take the data for the PG1115+080 system from Impey et al. (1998) (their 
Fig.l and Tables 1 and 2) who presented a near-infrared observation of the 
PG1115+080 system obtained with the Hubble Space Telescope (HST) NIC- 
MOS camera. The flux ratio of the close pair of images (Al and A2, see Fig. 
1) is approximately 0.67 and showed little variation with wavelength from 
the multiple wavelength observations by Impey et al. (1998) 0- Simple lens 
models require instead an A2/A1 flux ratio close to 1, because the images are 
symmetrically arranged near a fold caustic. There is no smooth lens model 
that can explain this anomalous flux: while each of such models can differ in 

2 Pooley et al. (2006), however, report that there has been some variation also in 
the optical. 
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complexity or in parameterization, the observed discrepancy in the flux ratio, 
compared with the expected universal relations for a cusp or fold singularity, 
suggests that it is an intrinsic difficulty for smooth lens models, not associated 
with a particular choice of the parameters (Yoo et al. 2005). 

Recently, Chiba et al. (2005) analyzed observations of the PG1 115+080 sys- 
tem done in the mid-infrared band and found a flux ratio A2/A1 of 0.93 for the 
close pair, which is virtually consistent with smooth lens models but clearly 
inconsistent with the optical fluxes. The observations indicated that the mea- 
sured mid-infrared flux originate from a hot dust torus around a QSO nucleus. 
Based on the size estimate of the dust torus, they placed limits on the mass of 
the substructure causing the optical flux anomaly! 3 1. For a substructure mod- 
eled as a SIS the subclumps should have a mass of at most 2.2 x 1O 4 M inside 
a radius of lOOpc to prevent anomalies in the mid-infrared band. However, it 
has to be pointed out that this latter result is based on several assumptions 
and few observations (Minezaki et al. 2004), so that the given value may be 
subject to substantial modifications. Indeed, if the size of the cooler dust torus 
causing the mid-infrared flux is underestimated then the above limit gets in- 
creased. Furthermore, Pooley et al. 2006 analyzed the system using recent 
X-ray observations, which show also a strong anomalous flux ratio. They do 
not exclude the microlensing hypothesis in order to explain the anomaly in the 
X-ray band, nonetheless they conclude that the optical emission region should 
be much larger (by a factor « 10 — 100) than predicted by a simple thin ac- 
cretion disk model. Within this model the source size should be R s ~ 10 15 cm 
(e.g. Wambsganss, Schneider, & Paczyhski 1990; Rauch & Blandford 1991; 
Wyithe et al. 2000). Therefore, if it is 10-100 bigger, R s « 0.01 - 0.1 pc, the 
effect of stellar microlensing could be ruled out (Metcalf 2005). 

An interesting quadruply imaged lens system is B1555+375 with a maximum 
separation of only 0.42 arcsec, which was discovered some years ago (Marlow et 
al. 1999). It has an anomalous flux ratio in the radio: the system was observed 
at 8,4 GHz at VLA and with MERLIN 5 GHz snapshot observations. There 
are only few observations in the optical and near-infrared band. Marlow et al. 
(1999) considered a model for B1555+375, which describes well the positions 
of the images but fails to reproduce accurately the flux ratio between the 
two images near the fold critical point (labeled by A and B, see Fig. 3). The 
observed ratio is about B/A ~ 0.57. This anomaly has also been discussed by 
Keeton et al. (2005) and Dobler & Keeton (2006). As the redshifts of lens and 
source have not yet been measured we will adopt the same values as used by 

3 They also considered for this system the microlensing hypothesis and estimated 
the time variability of the images flux in the mid-infrared band, which turned out 
to be rather long (more than a decade). This estimate is consistent with the fact 
that the optical flux ratio has remained unchanged over the past decade. It is thus 
clearly not yet possible to assess the microlensing hypothesis. 
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Image 



H I V 





mag 


mag 


mag 


Al 


15.75±0.02 


16.12 


16.90 


A2 


16.23±0.03 


16.51 


17.35 


B 


17.68±0.04 


18.08 


18.87 


C 


17.23±0.03 


17.58 


18.37 


Lens 


16.57±0.10 


18.40 





Table 1 " " 

Photometric data in 3 bands for the four images of PG1115+080, from Impey et al. 
(1998). 

Marlow et al. (1999): z t = 0.5 and z s = 1.5. 



3 Analytical treatment 



We briefly recall the general expressions for the gravitational lensing and refer, 
e.g., to the book by Schneider et al. (1992) (which we will denote afterwards 
with SEF) and the review by Kochanek (2004). The lens equation is 

= o- m , (i) 



where 0(9) is the source position and 9 the image position. a(9) is the deflec- 
tion angle, which depends on k(9) the dimensionless surface mass density or 
convergence in units of the critical surface mass density S C rit, defined as 



J crit 



c 2 Ds 
AttG D^Dls ' 



(2) 



where D$, D^, D^s are the angular diameter distances between observer and 
source, observer and lens, source and lens, respectively. 



3.1 Lens Mapping 



In the vicinity of an arbitrary point, the lens mapping can be described by its 
Jacobian matrix A: 

dfl ( da i (S)\ ( dH-m\ ,,, 
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Here we made use of the fact (see SEF), that the deflection angle can be 
expressed as the gradient of an effective two-dimensional scalar potential ip: 
a = Vgip, which carries information on the Newtonian potential of the lens. 
The magnification is defined as the ratio between the solid angles of the image 
and the source (since the surface brightness is conserved) and is given by the 
inverse of the determinant of the Jacobian A 



1 



detA' 



(4) 



The Laplacian of the effective potential ip is twice the convergence: 

ipn + ^22 = 2k = tr ipij. (5) 



With the definitions for the components of the external shear 7: 
H0) = l(^n ~ ^22) = 10) cos[2ip0)\ 



(6) 



and 



i20) = V12 = ^21 = 10) 



(7) 



(where the angle ip gives the direction of the shear) the Jacobian matrix can 
be written as 



A 




(1-k) 



' cos 2tp sin 2ip ^ 
^ sin 2ip — cos 2ip j 



(8) 



(9) 



where 7 = ylf + li- With eq.(8) the magnification can be expressed as a 
function of the convergence k and the shear 7 at the image point: 



H = (det^l)" 1 = 



1 



(1 - k) 2 - 7 2 ' 



(10) 



Locations at which det A = have formally infinite magnification are the 
critical curves in the lens plane. The corresponding locations in the source 
plane are the caustics. For spherically symmetric mass distributions the critical 
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curves are circles, whereas for elliptical lenses or spherically symmetric lenses 
with external shear, the caustics can have cusps and folds. 

Near a fold the lens equation can be reduced to a one-dimensional model and 
a Taylor expansion can be performed (see SEF, Kochanek 2004), for which we 
get 

^° = § (*-*o) + ^(0-0o) 2 , (11) 

i.e. 

(3 = 6{i-r)-\ro 2 ^-\re 2 (12) 

and inverse magnification 

p- 1 = (i - v") - i)'"e -> -i)"'e. (13) 

We choose the coordinate system such that there is a critical line at 9 = 
(i.e. 1 — ip" = 0) and the primes denote derivatives with respect to 9. These 
equations are easily solved and one finds that the two images are at 0± = 
±(—2P/ijj"') 1 ^ 2 with inverse magnifications /j,^. 1 = =p( — 2/5^ w ) 1//2 that are equal 
in magnitude but with opposite sign. Hence, if the assumptions for the Taylor 
expansion hold, the images merging at a fold should have identical fluxes. 
Using gravity to produce anomalous flux ratios requires terms in the potential 
with a length scale comparable to the separation of the images to significantly 
violate the rule that they should have similar fluxes. 

3.2 Perturbing the system 

Let's consider a general lens system configuration for which we know flux ratios 
and image positions and we assume to be able to reproduce with a smooth lens 
model, such as a singular isothermal ellipsoid (SIE), the main features of the 
lens, besides the anomalous flux ratio. Adding an external potential term in 
the lens equation and correspondingly in the Jacobian matrix such as induced 
by singular isothermal sphere (SIS) substructures perturb the system. Keeton 
(2001) analyzing the system B1422+231 (cusp case), could put limits on the 
sub clump mass range by considering the different effects they would induce 
on optical and radio bands. For the same system Bradac et al. (2002) suggests 
a way to estimate the minimum value for the convergence k in order to get 
agreement with the observed flux ratios. 
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Here we constrain the mass and the position of a substructure by considering 
its effects on the flux of the images. At each image position the perturbed 
Jacobian matrix can be written as 



1 - m - 71 



-72 



\ 



A = 



(14) 



-72 



1 



«i - 71 / 



where re x = (re + Are), 71 = (71 + A71), 72 = (72 + A72) and Are, A71 and A72 
are the perturbations induced by a substructure. 

If the substructure is modeled by a SIS, it is possible to express the shear 
components as a function of Are (this is true for models that have radial 
symmetry, Kormann et al. (1994)): A71 = AK,cos9 sis and A72 = AK,sin9 S i S , 
where 



(x S is, Usis) is the position of the substructure and (xp, yp) is the considered 
image position. R sis is the Einstein radius of the substructure, which depends 
on its mass and the distances (or redshifts) to the lens and the source, the 
latter ones being known quantities. The 9 S i S is given through the relation 
tg0 S i S = (x S i S — xp)/(y S i S — yp) and it is the angle between the SIS and the 
considered image position (xp,yp). We first consider a model with one ad- 
ditional substructure located at the same distance as the lens. Its mass and 
position have to be determined such that the substructure does not signifi- 
cantly (within the observational errors) modify the positions of all the images 
as well as the fluxes of the images which lie far from the two ones near the 
fold critical point. These requirements clearly put constraints on the mass and 
position of the substructure. 

For the determination of the magnification of an image the additional terms 
due to the substructure depend only on the position (x sis ,y sis ) and the mass 
(Einstein radius R S i S ) of the subclump, thus we have 3 unknown quantities 
(see Appendix). We consider only three images, thus getting a system with 
three equations for three unknown quantities, and assume that the 4th image 
is far enough such as not to be perturbed by the substructure. We then verify 
a posteriori that the found solution satisfies this latter assumption within the 
measurement errors. It turns out indeed to be the the subclump is 

located far from the 4th image, which is chosen as being the most distant 
one from the two near the fold. The system of equations is non-linear, so 
that the solution is not unique (see Appendix A). However, some solutions 
have to be discarded as being not physical (imaginary values or a negative 



Are 



■sis 



(15) 




9 



Einstein radius). All acceptable solutions are taken as input parameters for 
the simulations, as will be discussed in the next Section (see also Fig. 1). 

Note that the substructure could produce further multiple images of the orig- 
inal one. In our case, where we model the substructures as SIS, necessary and 
sufficient condition for multiple images formation is that the Einstein radius 
(of the subclump) 0Esub has to be greater than half of 9i, the distance of the 
image from the subclump (0Esub > (1/2)$/) (Narayan & Schneider 1990). In 
Tables 4 and 7 we give the positions of the substructures and of the images for 
the two considered lens systems as obtained from the simulations as discussed 
in the next section. From these data one easily verifies that the Einstein radius 
of the substructures, as given in Tables 2 and 6, do not satisfy the above con- 
dition, thus ruling out the formation of further images. As noticed by Keeton 
(2003), for SIS subclumps positive-parity images get always brighter, whereas 
negative-parity images get fainter. In our case for PG1115+080 (B1555+375) 
Al (A) is the positive-parity image and A2 (B) the negative-parity one. 



4 Numerical simulations 



In this section we present our simulations and results. We use the gravlens 
code developed by Keeton (20010, modeling the main galaxy acting as lens 
(in both cases) by a SIE and then by adding an external shear term (which we 
will in the following denote by SIE^) and/or a SIS term to take into account 
the influence of the group in which the galaxy is embedded. Moreover, we add 
one or two substructures to take into account the effects on small scales. For 
PG1115+080 we use data in H band (see Table 1) taken from Impey et al. 
(1998), and for B1555+375 the data in the 5 GHz radio band (see Table 6) from 
Marlow et al. (1999). We allowed a conservative 1 a error in the relative x- and 
y- positions of the image components (corresponding to an error of at most 5 
mas), and 1 a error on the values of the fluxes (corresponding to a variation 
by 20%). In each model we have different parameters and constraints, and the 
goodness of the fit is given by the x 2 value, evaluated on the image plane 
{Ximg)f an d is a sum °f different contributions: image positions and fluxes, and 
main galaxy position^- 



4 The software is available via the web site: |http: / / cfa-www .harvard. edu /castles| 

5 There is an alternate way to define the \ 2 that is evaluated in the source plane 
(xLc) ( e -S-> Kayser et al. 1990), which is an approximate version ofxf mg ' when using 
the minimization within this approximation the formation of additional images is 
not excluded, maybe yielding to a not realistic model. However, the approximation 
inherent Xsro should properly be used only if a good model is already known, not 
in an initial search for a good model (Keeton 2001). 
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Fig. 1. PG1115+080: images, source and galaxy (G) positions are shown assuming a 
SIE 7 +SIS model. The position of the substructure for both models is also given: note 
that SI is closer to the Al image, while S2 to A2. The critical curve and the caustic 
are for the SIE 7 model without the modifications induced by the substructure. The 
analytic solutions as discussed in the Appendix are also shown. The solutions labeled 
as 3 and 4 are so close that on the figure they coincide. 

4.1 PG1115+080 



We model the main lens galaxy as a SIE and take into account the contribution 
due to the group of which the galaxy is part of either by adding a SIS group 
component (Keeton 2003, Chen et al. 2003) or an external shear term. We 
consider both models to which we add one substructure described as a small 
SIS. 

The results for PG1115+080, modeled as a SIE and an external shear or a 
SIS, are given in Tables 2, 3 and 4 ( in the following we indicate with M the 
mass inside the Einstein radius) and in Fig. 1. In both cases we find good 
agreement with previous results (Impey et al. 1998, Chiba 2002). In the case 
SIE~f we find for the external shear 7=0.11 and for its direction an angle of 
= 56°, which agree well with the results of Chiba (2002). 

In a further step we add one substructure, using as starting parameters for 
its position and mass the analytic results as determined following the method 
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outlined in the previous section 0. In the model SIE 7 + SIS we have 11 
parameters, i.e. for the main galaxy : Einstein radius (i.e. mass), ellipticity 
e and orientation PA, shear 7 and its direction 0; for the substructure: the 
Einstein radius, position (corresponding to 2 parameters) and, moreover, the 
source position (2 parameters) and its flux. The observational constraints are 
12, namely the 4x2 image positions and 4 fluxes. 

In the model SIE+SIS 9rcmp +SIS, in which the group is modeled by a SIS, the 
substructure (denoted in Fig.l by S2) is close to A2, whereas for the previous 
model it is close to Al (denoted as SI in Fig. 1). By adding one substructure 
the value of the anomalous flux ratio improves substantially for both models: 
getting lowered from 0.91 to 0.69 for the first model and to 0.66 for the second 
model, with xj ot =\.?> and « 0.4, respectively, which is quite good. 

The results from the simulations agree, as expected, quite well with the an- 
alytical ones. Since the mass of the substructures is very small as compared 
to the mass of the lens galaxy, the approximation used in the analytical ap- 
proach to neglect the influence induced on the positions of the images by the 
substructures is quite well fulfilled. We checked this, indeed, on the results 
obtained from the numerical simulations. 

Since the positions and magnifications of the images are only known within 
a certain accuracy, we computed the corresponding la and 2a ranges for the 
value of the mass of the substructure. Starting from our best model we consider 
two approaches. In the first we let all parameters (i.e. the main galaxy ones 
and the position of the substructure) vary, whereas in the second one we keep 
the main galaxy parameters fixed at the values given by the best fit model and 
let the remaining parameters vary. The high values for the total x 2 are due to 
the bad galaxy position fit (in the first case) and to the bad image position fit 
(in the second case). In Fig. 2 we report x 2 as a function of the substructure 
Einstein radius for the SIE+SIS 9r . oup -|-SIS model for both cases mentioned 
above. We do not consider larger values for the Einstein radius as this would 
correspond to masses for the substructure of order M{< Re) ~ 10 9 — 1O 1O M , 
too big to be realistic and for which one would see effects on the image position 
or even image splitting. On the other hand we can also exclude Einstein radius 
that are too small. In fact, if it is true that with a standard accretion disk 
model we get a source size radius R s ~ 10 15 cm (Chiba et al 2005) and that 
for the PG1115+080 system the real source size should be about 10-100 times 

6 Moreover, since we don't know a priori the source flux, we take the value we get 
from the SIE^ model and, considering the observational fluxes with respect to the 
Al image, we use them as starting values in the analytical system. 

7 By further adding a second substructure to the first model, we find even a lower 
value for the flux ratio, however, Xtot an d Xfi ux increase, which indicates that this 
model does not correspond to a global minimum, but rather to a local minimum of 
the Xtot surface. 
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Parameter 


SIE 7 


SIE+SISg rcm p 


SIE 7 +SIS 


SIE-hSIS^f^p-hSIS 


RE,gal 


1.03 


1.14 


1.12 


1.03 


M gal (M Q ) 


1.23 x 10 11 


1.05 x 10 11 


1.24 x 10 11 


1.05 x 10 11 


Re , group 


— 


2.30 


— 


2.11 


M groU p (M ) 


- 


4.10 x 10 11 


- 


4.0 x 10 11 


RE,subl 


— 


— 


0.033 


0.001 


M subl {M Q ) 


— 


— 


1.00 x 10 8 


1.00 x 10 5 


<jS IE (kms~ l ) 


232.3 


245.4 


243.5 


232.3 




— 


— 


39.2 


6.8 


A2/A1 


0.91 


0.95 


0.69 


0.66 


7 


0.11 




0.11 




4> 


56° 




56° 




x 2 


77 


3.9 


1.30 


0.08 



Table 2 



PG1115+080: Parameters for different models, with and without substructure. 7 
and 4> are the values of the external shear and its direction. Re is expressed in arcsec 
(1 arcsec = 3.19 kpc h~ l ). The ID velocity dispersions for the main component and 
substructures are also reported. By adding a substructure the agreement between 
predicted and observed fluxes increases substantially. 



Model 


Image 




K 




7 




A2/A1 


SIE 7 


Al 


0. 


.498 





.421 


13.35 


0.91 




A2 


0. 


535 





.545 


-12.17 




SIE+SISg rcm p 


Al 


0. 


.534 





.411 


20.21 


0.95 




A2 


0. 


.551 





.504 


-19.31 




SIE 7 +SIS 


Al 


0. 


.554 





.372 


16.54 


0.69 




A2 


0. 


.561 


0. 


.531 


-11.54 




SIE+SISg roU p+SIS 


Al 


0. 


.531 





.410 


19.46 


0.66 




A2 


0. 


.565 





,517 


-12.83 





Table 3 

PG1115+080: Values of shear, convergence and amplification for Al and A2 images 
from simulations for the considered models. 
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Ubject 


X 


y 


„ DA 

e JrA 




(arcsec) 


(arcsec) 




Galaxy 


-0.381 


-1.345 


n 1 a a 
(J. 14 — 84 


Al 


-1.328 


-2.037 




A2 


-1.478 


-1.576 




Subl 


-1.33 


-2.20 




Sub2 


-1.52 


-1.57 





Table 4 

PG1115+080: positions of the lens galaxy center, the close pair Al and A2, as well 
as the substructure with respect to the C image (see Fig.l). Also ellipticity e and 
orientation PA of the semi major axis with respect to x-axis (as measured from 
East to North and centered in the C image) are given. The distances between the 
substructure and the images Al and A2 are bigger than twice their corresponding 
Einstein radius. Thus no further images will be formed. 

bigger (Pooley et al 2006), we can estimate roughly the limit of the Einstein 
radius for which the effects on the images become negligible. For a stellar 
Re ~ 10 -5 arcsecs, that corresponds to 0.03 pc (on the lens plane), there 
would be no (or little) effect on an image of a source with Rs ~ 0.1 pc. The 
minimum value of the curve corresponds to an Einstein radius ~ 0.001 arcsec. 
Anyway, the curve within the 2a range is rather constrained (in both cases), 
leading thus in practice to a small degeneracy, with a rather narrow range 
for acceptable values of the substructure mass. The 2a range is within an 
Einstein radius of 0.0005 and ~ 0.005 arcsec corresponding to ~ 2.5 x 10 4 M 
and ~ 2.5 x 10 6 M . 

The 2a range for the SIE^ + SIS model is somewhat larger. Considering for 
example the first case, the Einstein radius for the substructure could be as 
high as 0.1 arcsecs, leading to quite a big mass (~ 1O 9 M ). 



4.2 B 1555+375 

B1555+375 is another lens system for which the agreement between observa- 
tions and model can be improved. Already by adding an external shear one 
can quite substantially improve the flux fitting, assuming for the lens galaxy a 
SIE model. Anyway, Dobler & Keeton (2006) consider this model unphysical 
(since ellipticity and shear turn out to be almost perpendicular) and discuss 
other models, adding substructures and giving lower limits on their masses. 
They find an acceptable model using two substructures in front of B and C 
images, respectively. Adding one substructure to the model SIE 1 does not 
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Fig. 2. la and 2a limits for the Einstein radius for the model SIE+SIS group +SIS: 
a): letting all parameters fixed but the position of the subclump; b): letting all 
parameters vary. 

improve the fit furtheipH On the other hand, adding two substructures (mod- 
eled as SIS with comparable masses and located near the pair of images due 
to the fold) to the simple SIE model, we get a value for the flux ratio B/A of 
0.59, which is only slightly higher then the observed one of 0.57 The results 
are reported in Tables 6, 7 and 8 and in Fig. 3. For the model of a SIE 1 we 
get x 2 — 1-29 and xjiux = 0.35. For the solution with two substructures we 
find a x 2 ~ 5, with N do f = 5. 

As we pointed out above, in the radio band up to now the only successful 



We tried, as in the previous case, also analytically to estimate mass and position 
for a single SIS added to a SIE model, but the system of equations does not have 
in this case real solutions. 
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Image 


X 


y 


S 5 




(arcsec) 


(arcsec) 


(mJy) 


A 


0.0±0.005 


0.0±0.005 


17.0 


B 


-0.073±0.005 


0.048±0.005 


9.7 


C 


-0.411±0.005 


-0.028±0.005 


8.3 


D 


-0.162±0.005 


-0.368±0.005 


1.3 



Table 5 ' ' ' 

B1555+375: Positions and photometric data of the 4 images as given by CLASS 
(from Marlow et al. 1999). 

explanation for the flux anomalies is to consider substructures. Moreover, we 
notice that the SIE model does neither fit very well the fluxes of the other 
images besides the close pair ones. 

As mentioned in the SIE^ model we get a rather high value for the ellipticity 
(0.85) and for the external shear (0.23): such a strong shear can be induced 
by a group of galaxies located around the main lens. One has also to con- 
sider possible effects due to groups of galaxies which lie on the line-of-sight, 
both in the foreground and in the background (Williams et al. 2006). Other 
systems (like B1608+656 or HST 12531-2914) show such a high value for the 
external shear (see Witt & Mao 1997). In particular, B0128+437 seems to be 
quite similar to B1555+375 (see Philipps et al. 1999). Anyway, there are still 
not enough observations neither about the main galaxy nor about its envi- 
ronment, so that it is not possible to choose among the different solutions. 
However, our solution with two substructures involves a SIE lens model with 
an angular structure which is in agreement with previous works (Marlow et 
al. 1999 and Keeton, Gaudi & Petters 2005). Also for this model we computed 
the confidence intervals for the various parameters. As an example in Fig. 4 
we show the contour ellipses for the confidence intervals of the Einstein radii 
of the substructures as obtained fixing the main galaxy parameters. 



5 Discussion 

We have analyzed two lens systems, PG11 15+080 and B1555+375, which show 
an anomalous flux ratio for the two images near the critical curve, due to a 
fold configuration. These systems cannot be modeled using only smooth lens 
models like SIE, although they fit well all the positions of the images. We added 
one or two substructures, taking as starting parameters for our numerical 
simulations the ones obtained by an approximated analytical treatment. 

In this way we reproduce well, in addition to the positions, also the fluxes 
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Parameter 


SIE 


SIE 7 


SIE+2SIS 


RE,gal 


0.22 


0.165 


0.21 


Mg a l (M ) 


0.50 


0.45 


0.50 x 10 10 


RE,subl 


— 


— 


0.009 


RE,sub2 


— 


— 


0.012 


M subl (M Q ) 


— 


— 


8.1 x 10 6 


M sub2 (M Q ) 






1.4 x 10 7 


of ^(ferns' -1 ) 


170.1 


147.3 


165.6 


a^ bl (ferns -1 ) 


— 


— 


23.3 


a s v ab2 (kms- 1 ) 


— 


— 


24.7 


B/A 


0.93 


0.61 


0.59 


e 


0.53 


0.85 


0.53 


PA 


-2.42° 


-6.39° 


0.41° 


7 




0.23 








-78° 




x 2 


13.9 


1.9 


5.2 



Table 6 

B1555+375: Results from the simulations for two models without substructures 
and one model with two substructures. (Notice that for the latter model the \ 2 
is higher, see text). The Einstein radii are expressed in arcsec. The system is well 
fitted already by adding external shear. 



Object 


X 

(arcsec) 


y 

(arcsec) 


e 


PA 


Lens 


-0.162 


-0.246 


0.53 


0.75° 


A 


0.0 


0.0 






B 


-0.075 


0.043 






Subl 


-0.060 


0.094 






Sub2 


0.101 


0.001 







Table 7 

B1555+375: Parameters of the lens model and of the added substructures, e and 
PA are ellipticity and orientation of the semi major axis with respect to x-axis (as 
measured from East to North and centered in the A image). 
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Fig. 3. B1555+375: images, source and galaxy (G) positions are shown assuming a 
SIE+2SIS model. The positions of the two substructures are also given. The critical 
curve and the caustic are for the SIE model alone without the modifications induced 
by the two substructures. 

of the pair of images near the critical curve. In the PG1115+080 case we 
get a ratio of 0.69 as compared to the observed one of 0.65 if we consider a 
model SIE 7 + SIS or 0.66 for a model SIE+SIS 5roup +SIS and for B1555+375 
a ratio of 0.59 instead of 0.57 with a SIE+2SIS model. We point out that for 
PG1115+080 when we consider the model SIE+SIS 9rcmp + SIS we find that the 
substructure needed to explain the anomaly in the optical band lies close to 
the A2 image and has a mass ~ 10 5 M©. On the other hand, the model SIE 7 
10 8 M Q close to the Al image that could in principle affect 
also larger A bands, for which no anomaly has been reported yet (Chiba et 
al 2005). Therefore, this latter model is certainly less plausible if not already 
excluded. Clearly, new observations are still needed to better constrain the 
models. 

In the case of B1555+375 we also explored simpler models (i.e. SIE, SIE plus 
simple external shear, or SIE + SIS) but we did not find any acceptable 
solution. We note that our best model is similar to the one found by Dobler & 
Keeton 2006: in their analysis, they concluded that the system B1555+375 has 
two anomalous images. Their approach is different, since they try to identify 
anomalies by relaxing the flux constraints on an image and fit all the others, 
image positions and fluxes. When a good model is found, they than conclude 
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Fig. 4. B1555+375: Contour plot for the 1,2 and 3a confidence intervals for the 
Einstein radii (in arcsec) of the substructures. 

that the unconstrained image is anomalous and give a lower mass limit on 
the sub clump likely responsible for the anomalous flux by optimizing \ 2 as 
a function of the subhalo Einstein radius and the source size. Nevertheless, 
they also need two substructures close to two different images (either in front 
of images A and C or close to images B and C) to fit the fluxes. Despite the 
fact that one of their models has N^of = 0, the values of the subhalo masses 
are similar to ours (e.g. ~ 10 5_6 M Q within the Einstein radius). Finally, we 
would like to point out that the current observations towards B1555+375 seem 
not to reveal any nearby group to which the lens system belongs, so a SIE 
model plus external shear might be inadequate to describe the effect due to 
substructures, which are rather better modeled by a SIS. 

In both cases the best models are the ones with additional substructures. 
Previously, the best models gave values for the anomalous flux ratio of 0.91 
for PG1115+080 and 0.93 for B1555+375, respectively. The improvement 
achieved by just adding substructures is remarkable. 

Our approach, although leading to better values for the anomalous flux ratios, 
is not completely new and indeed, as mentioned in the introduction previous 
works attempted to solve the anomalous flux ratios problem by adding per- 
turbations until a good model was found. In particular, Keeton 2003, Inoue 
& Chiba 2005 studied the effect of substructures, modeled as SIS, in a con- 
vergence and shear field describing the main lens galaxy (or galaxies along 
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the line of sight, e.g. Keeton 2003). Mao & Schneider 1998, Metcalf & Madau 
2001, Keeton 2001 analysed the problem by directly adding subclumps, in 
form of point masses, SIS, or plane- wave perturbations, in order to explain 
the observed anomalies. Chen 2003, Metcalf 2005, Keeton et al. 2005, Maccio 
& Miranda 2006 (and references therein) studied the problem by using statis- 
tical approaches and numerical simulations. 

The masses of the substructures are in the range w 10 5 - 10 8 M Q , and distant 
enough from the images not to induce the formation of additional ones. Given 
the mass range the perturbers could be globular clusters or small satellite 
galaxies, but also CDM dark substructures. In order to compare the masses 
we obtained with dark substructures as predicted by CDM, we notice that our 
values correspond to the mass enclosed within the Einstein radius, while CDM 
subhalo masses are defined as the mass enclosed within the tidal radius, which 
is usually much larger than R E . However, without going too much into details, 
it is still possible to get a rough estimate of the mass within the tidal radius 
by calculating the ID velocity dispersion of the subhalos from their Einstein 
radius and then by considering the M oc v f (with the approximation v c = 
&v/V3) relation found in numerical simulations (Bullock et al. 2001, Diemand 
et al 2004). Althought current simulations achieve a lower limit for CDM 
subhalo masses of M w 1O 1O_11 M , we assume this relation to hold for lower 
values of the velocity dispersion and of the masses. Clearly, the extrapolation 
to lower velocities could suffer from resolution effects and numerical noise. In 
the system PG1 15+080 for the SIE 1 + SIS model the substructure has a 
a v w 39km/ s, which leads to a mass of ~ 7 x 1O 8 M (the mass within the 
Einstein radius, as given in Table 2, is instead 1 x 1O 8 M ), while in the model 
SIE + SIS group + SIS the substructure has a v ~ 6.8km/ s corresponding to ~ 
4 x 1O 6 M (instead of 1 x 10 5 M Q ). In the system B1555+375 the substructures 
have a v i ps 23 km/s and cr v2 ~ 25 km/s respectively, which leads for both a 
corresponding mass of about 1O 8 M (whereas, as given in Table 6, the mass 
within the Einstein radius is ps 1O 7 M ). We roughly find that the mass within 
the tidal radius is about 10 times bigger than the one within the Einstein 
radius. 

To get a rough estimate of the number of substructures expected to lie close 
to the images we follow the work of Diemand et al. (2004). They compute 
the two dimensional radial number density of subhalos inside a galaxy virial 
radius, from which it is then possible to get an estimate of the number of 
substructures inside a small area surrounding an image in a lens system. The 
number of subhalos with a mass greater than m inside an area A at a distance 
r from the center of the galaxy is given by (Maccio & Miranda 2006) 

N A (> m, r) = " ;/ 2 m V ; , 16 
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where {N Tv (> mo)) is the average number of subhalos with m > itlq inside the 
virial radius r v of the galaxy and N(r) describes the radial dependence of the 
number of substructures. The cumulative mass function of subhalos within 
the virial radius of an halo scales as oc mr 1 . 

As an example for PG1115+080, we consider m > 10 5 M Q , (mo ~ 10 7 ) and 
a distance of the images from the center r rs 1.5 arcsec (N(r) m 2 — 6, 
(N rv (> m )) ~ 166,r vir = 268 kpc, see simulations GO and Gl in table 1 of 
Diemand et al. 2004 and Maccio & Miranda 2006). Typically, a substructure 
is expected to lie at a distance from an image comparable to the separation 
between the close pair, which is about ~ 0.5 arcsec. We, therefore, consider an 
area corresponding to a small disc with a radius about twice the separation 
distance, thus ~ 1.0 arcsec (corresponding at Z\ = 0.31 to Am 7r(4.5) 2 /cpc 2 ). 
With these assumptions we expect 10 to 30 substructures in the considered 
area. If instead we require m > 1O 6 M then we find about 1 to 3 substructures, 
and 0.1 to 0.3 for m > 1O 7 M . From these considerations we see that the 
expected number of substructures within the CDM model seems, in particular, 
to be lower than required when considering the model with a 10 8 M Q mass 
subclump or even bigger, whereas for models with a 1O 6 M mass subclump 
it might be in better agreement. 

For B1555+375 the above considerations apply as well: at z% = 0.5, 1.0 arcsec= 
6.114 kpc, thus the area around the image pair is Am n (6. 114) 2 kpc 2 . By using 
the same values as in the previous case the expected number of CDM subhalos 
bigger then 10 8 being close to the image pair is between 0.02 and ~ 0.1. 
Clearly, given the rough approximation used, some of them depending on 
extrapolations of numerical simulations of limited resolution, it is not possible 
to draw firm conclusions. 

It is obviously also not possible to distinguish between different models, more- 
over it could be that some of the substructures are actually located along 
the line of sight rather than being in the surroundings of the lens galaxy. To 
this respect it is interesting to notice that with future ALMA observations 
one could solve this latter problem as pointed out by Inoue & Chiba (2005). 
They proposed a method to realize a 3D mapping of CDM substructures in 
extragalactic halos, based on astrometric shift measurements (at submillime- 
ter wavelengths) of perturbed multiple images with respect to unperturbed 
images, with which it should be possible to break the degeneracy between the 
subhalo mass and a position along the line of sight to the image. 

Also other explanations of the flux anomaly in multiple QSOlens systems have 
been considered in the literature, however the best solution, seems to be the 
presence of substructures in the halo of the lens galaxy. 

For the second case (B1555+375) there are two acceptable solutions: with and 
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without any substructure. Even if a high shear value seems unplausible, we 
cannot yet rule out this possibility. A recent work (Williams et al. 2006) shows 
the importance of galaxy groups along the line of sight that can significantly 
impact the lens model. High-resolution VLA radio observations could help to 
constrain the lens model further. On the other hand, starting from the simple 
SIE model, we find a real good fit for the image positions and fluxes, if we add 
two substructures, still keeping acceptable values for the main lens parameters 
(in agreement with Marlow et al. 1999). We note also that for the B1555+375 
system we assumed the redshifts of the source and the lens to be known. Even 
if these values might not be exact the uncertainties do not change the result 
significantly. For instance, if we let the redshifts (z s = 1.5 and zi = 0.5) vary 
by ±0.3 and use a source redshift of z s = 1.8 with z\ = 0.2 we get a mass for 
the first substructure of 3.3 x 1O 6 M (< Re)- By considering instead z s = 1.2 
and zi = 0.8 the mass is 2.7 x 1O 7 M (< R E ). All other combinations will give 
a mass value within this range, and similarly for the second substructure. 

Finally, we observe that the PG1115+080 system is radio-quiet, so that the 
microlensing hypothesis can not be ruled out (Pooley et al. 2006). However, 
given the different observations at the various wavelengths it might also be 
possible that both microlensing and millilensing are at work. A more accu- 
rate analysis about the source size could cast some light in constraining the 
substructure size. Since there are discussions about it, more high resolution 
observations are needed to definitely rule out the millilensing or microlensing 
hypothesis. On the other hand, for the B1555+375 system anomalies are evi- 
dent in the available radio data, for which the most plausible explanation are 
CDM substructures in galactic halos. 
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A Analytic estimates for convergence and shear due a substructure 

We briefly present here the formalism used for the analytical approximation 
of the total convergence K tot and the total shear / -f t ot in the presence of one 
perturber, located in the main lens plane. At each point of the lens plane we 
can evaluate the total amplification using the quantities: 
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Ke 


X 


y 




arcsec 


arcsec 


arcsec 


1 


U.UoU 


1 1 on 

-i.iyy 


-Z.Ob 


2 


0.035 


-1.825 


-2.352 


3 


0.049 


-1.363 


-2.421 


4 


0.050 


-1.76 


-2.332 



Table A.l 

Analytic solutions for PG1115+080 as discussed in the text. 



^tot — ^sie 4" Aft, 
Tltot = 7lsie + A71, 
72tot = l2sie + A7 2 , 

A*" 1 = (1 - Ktotf ~ (lltotf + (l2tot) 



(A.l) 
(A.2) 
(A.3) 

(A.4) 



n sie is the convergence due to the main lens and An is due to the perturber 
and similarly for the shear ^ sie and A7 (see Sec. 3.2). Dealing with a SIE 
model allows us to write (see Kormann et al. 1994): 



Re 



Kg, 



2q 2 
l+q- 



■){x sie -X P ) 2 + {Vsie-ypf 



(A.5) 



Where q is the axis ratio of the elliptical model used for the galaxy act- 
ing as lens. The values for K S i e , 7i S j e and r )2sie are taken, for instance in the 
PG1115+080 case, from the SIE^ model (which corresponds to the zero order 
approximation), computed in Al and A2 (see Table 3) as well as the magni- 
fication factors with respect to the image Al (and keeping the source flux as 
obtained from SIE^ model). In this way the system 

HAL = (1 - KtotAlY - [{lltotAl) 2 + {l2totAlf] 

< fl A \ = (1 - K totA2 ) 2 - [(TUtotAi) 2 + iri2totA2?] ( A -6) 
^B 1 = (! - K totB? ~ [{lltotBf + (l2totB) 2 } 

v 



has only three unknown quantities, namely R S i S and the perturber position 
given by (x sis , y S i S )- Since the system is non-linear, we get different sets of 
solutions, some of which turn out to be unphysical. The allowed solutions are 
near the close pair (see Table Al). We marked their positions on Figure 1. We 
used then these solutions as input parameters for the numerical simulation. It 
turns out that all converge to the same model SIE 1 + SIS discussed in Sect. 
4.1. 
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